W wyniku generacji układu równań dostaliśmy macierz zapisaną w postaci produktu Kroneckera dwóch macierzy pięcio-przekątniowych, oraz prawą stronę policzoną z mozołem dla poszczególnych B-spline'ów testujących. Musimy teraz rozwiązać uzyskany w ten sposób układ równań.
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{20} & \frac{13}{120} & \frac{1}{120} \\ \end{bmatrix} \otimes \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ u_{3,1} \\ \vdots \\ u_{k,l} \\ \vdots \\ u_{N_{x-2},N_y} \\ u_{N_{x-1},N_y} \\ u_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} \int BITMAP(x,y) B^x_1(x)*B^y_1(y) dxdy \\ \int BITMAP(x,y) B^x_2(x)*B^y_1(y) dxdy \\ \int BITMAP(x,y) B^x_3(x)*B^y_1(y) dxdy \\\vdots \\ \int BITMAP(x,y) B^x_k(x)*B^y_l(y) dxdy \\\vdots \\\int BITMAP(x,y) B^x_{N_x-2}(x)*B^y_{N_y}(y) dxdy \\\int BITMAP(x,y) B^x_{N_x-1}(x)*B^y_{N_y}(y) dxdy \\\int BITMAP(x,y) B^x_{N_x}(x)*B^y_{N_y}(y) dxdy \\ \end{bmatrix} \)
Rozwiązanie układu równań, w którym macierz ma strukturę produktu Kroneckera jest możliwe w bardzo krótkim czasie. Co to znaczy w bardzo krótkim czasie? Koszt obliczeniowy wyraża się liczbą operacji takich jak mnożenie czy dodawanie liczb, koniecznych do rozwiązania układu równań. W przypadku układu równań, w którym macierz ma strukturę produktu Kroneckera możliwe jest rozwiązanie układu równań za pomocą algorytmu, w którym liczba operacji wynosi \( const*N \) gdzie \( N \) to jest liczba niewiadomych (liczba współczynników aproksymacji bitmapy na siatce, wyrażona dokładnie poprzez \( N=N_x*N_y \), gdzie \( N_x,N_y \) to rozmiary siatki, natomiast \( const \) oznacza pewną stałą liczbę.
Algorytm stosowany w tym przypadku nazywany jest algorytmem solwera zmiennokierunkowego. Rozważamy dwa etapy procesu rozwiązania. Pierwszy etap polega na wzięciu pierwszej z podmacierzy, oraz poukładaniu wektora niewiadomych i wektora prawych stron w wiele podwektorów, po jednym wektorze dla każdej kolumny elementów siatki obliczeniowej. Zostało to zilustrowane na Rys. 1, oraz we wzorze poniżej.
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} w_{1,1} & w_{1,2} & \cdots & w_{1,N_{y-1}} & w_{1,N_y} \\ w_{2,1} & w_{2,2} & \cdots & w_{2,N_{y-1}} & w_{2,N_y} \\ w_{3,1} & w_{3,2} & \cdots & w_{3,N_{y-1}} & w_{3,N_y} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ w_{N_{x-2},1} & w_{N_{x-2},2} & \cdots & w_{N_{x-2},N_{y-1}} & w_{N_{x-2},N_y} \\ w_{N_{x-1},1} & w_{N_{x-1},2} & \cdots & w_{N_{x-1},N_{y-1}} & w_{N_{x-1},N_y} \\ w_{N_x,1} & w_{N_x,2} & \cdots & w_{N_x,N_{y-1}} & w_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} \int BITMAP(x,y) B^x_1(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_1(x)*B^y_{N_y}(y)dxdy \\ \int BITMAP(x,y) B^x_2(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_2(x)*B^y_{N_y}(y)dxdy \\ \int BITMAP(x,y) B^x_3(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_3(x)*B^y_{N_y}(y)dxdy \\ \vdots \\ \int BITMAP(x,y) B^x_{N_x-2}(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_{N_x-2}(x)*B^y_{N_y}(y)dxdy \\ \int BITMAP(x,y) B^x_{N_x-1}(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_{N_x-1}(x)*B^y_{N_y}(y) dxdy \\ \int BITMAP(x,y) B^x_{N_x}(x)*B^y_1(y) dxdy & \cdots & \int BITMAP(x,y) B^x_{N_x}(x)*B^y_{N_y}(y)dxdy \\ \end{bmatrix} \)
Wprowadziliśmy tutaj pomocnicze niewiadome, \( w* \) które służą do rozwiązania pierwszego układu równań. Oryginalne niewiadome \( u* \) dostaniemy po rozwiązaniu drugiego układu równań, w którym prawą stronę stanowić będą niewiadome pomocnicze \( w* \). Dostaliśmy więc układ równań z macierzą pięcio-przekątniową, o wielu prawych stronach. Każdy podwektor, każda prawa strona, odpowiada jednej kolumnie na siatce elementów, ma więc ustaloną współrzędną \( y \), oraz współrzędną \( x \) zmieniającą się od 1 do \( N_x \).
Podobnie uporządkowane są niewiadome \( u* \), w których to wierszami zmieniają się drugie indeksy, na przykład \( w_{1,1}, w_{1,2}, ..., w_{1,N_y} \) natomiast w kolumnach zmieniają się pierwsze indeksy.
Rozwiązujemy ten układ równań (jak to zrobić o tym za chwilę), dostajemy rozwiązania \( w* \) i przechodzimy do drugiego kroku algorytmu solwera zmiennokierunkowego.

Drugi etap polega więc na wzięciu analogicznej drugiej macierzy produktu Kroneckera, wzięciu rozwiązań \( w* \) z pierwszego układu równań, poukładaniu ich zgodnie z wierszami elementów z siatki (porównaj Rys. 1, poukładaniu w podobny sposób poszukiwanych niewiadomych \( u* \) i na rozwiązaniu otrzymanego układu równań o wielu prawych stronach.
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\\cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} u_{1,1} & u_{2,1} & \cdots & u_{N_{x-1},1} & u_{N_x,1} \\ u_{1,2} & u_{2,2} & \cdots & u_{N_{x-1},2} & u_{N_x,2} \\ u_{1,3} & u_{2,3} & \cdots & u_{N_{x-1},3} & u_{N_x,3} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ u_{1,N_{y-2}} & u_{2,N_{y-2}} & \cdots & u_{N_{x-1},N_{y-2}} & u_{N_{x},N_{y-2}} \\ u_{1,N_{y-1}} & u_{2,N_{y-1}} & \cdots & u_{N_{x-1},N_{y-1}} & u_{N_{x},N_{y-1}} \\ u_{1,N_y} & u_{2,N_y} & \cdots & u_{N_{x-1},N_{y}} & u_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} w_{1,1} & w_{2,1} & \cdots & w_{N_{x-1},1} & w_{N_x,1} \\ w_{1,2} & w_{2,2} & \cdots & w_{N_{x-1},2} & w_{N_x,2} \\ w_{1,3} & w_{2,3} & \cdots & w_{N_{x-1},3} & w_{N_x,3} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ w_{1,N_{y-2}} & w_{2,N_{y-2}} & \cdots & w_{N_{x-1},N_{y-2}} & w_{N_{x},N_{y-2}} \\ w_{1,N_{y-1}} & w_{2,N_{y-1}} & \cdots & w_{N_{x-1},N_{y-1}} & w_{N_{x},N_{y-1}} \\ w_{1,N_y} & w_{2,N_y} & \cdots & w_{N_{x-1},N_{y}} & w_{N_x,N_y} \\ \end{bmatrix} \)
W tym drugim układzie równań, każdy podwektor, każda prawa strona, odpowiada jednemu wierszowi na siatce elementów, ma więc ustaloną współrzędną \( x \), oraz współrzędną \( y \) zmieniającą się od 1 do \( N_y \). Podobnie uporządkowane są niewiadome \( u* \), w których to wierszami zmieniają się pierwsze indeksy, na przyklad \( w_{1,1}, w_{2,1}, ..., w_{N_x,1} \) natomiast w kolumnach zmieniają się drugie indeksy.
Każdy z tych dwóch układów równań o wielu prawych stronach rozwiążemy używając algorytmu eliminacji Gaussa dla macierzy pasmowej, który ma liniowy koszt obliczeniowy.